2. Análisis

Una vez realizado el preprocesado de los datos, realizamos un análisis de los datos obtenidos. Para ello, en primer lugar cargaremos las librerias que vamos a utilizar y los datos obtenidos del preprocesado.

library(readr)      
library(data.table) 
library(utils)     
library(sf)
library(dplyr)
library(purrr)
library(lubridate)
library(ggplot2)
library(corrplot)
library(factoextra)
library(leaflet)
library(leaflet.extras)
df <- read_csv("data_complete.csv", locale = locale(encoding = "UTF-8"))

1. Preparación de los datos

Para un análisis de la distribución de las variables, en primer lugar prescindiremos de las variables geométricas dentro del dataframe:

# Descartamos los valores geométricos
df <- df %>% select(-geometry, -geom)

A continuación veremos como es la estructura interna de los datos importados:

str(df)
tibble [20,065 × 37] (S3: tbl_df/tbl/data.frame)
 $ ...1                  : num [1:20065] 1 2 3 4 5 6 7 8 9 10 ...
 $ FID                   : num [1:20065] 6188 6572 6589 6606 6695 ...
 $ admin_region          : num [1:20065] 6 4 4 3 2 2 6 6 1 3 ...
 $ admin_district        : num [1:20065] 2 2 8 6 9 9 1 2 4 1 ...
 $ municipality          : num [1:20065] 603 405 418 318 221 221 601 603 111 306 ...
 $ year                  : num [1:20065] 2016 2016 2016 2016 2016 ...
 $ month                 : num [1:20065] 4 2 3 3 1 3 2 3 3 4 ...
 $ hour                  : num [1:20065] 15 15 14 13 8 7 7 15 15 16 ...
 $ day_week              : num [1:20065] 2 4 6 7 6 5 5 3 7 3 ...
 $ category              : num [1:20065] 3 3 3 3 3 3 3 3 3 3 ...
 $ kind_accident         : num [1:20065] 5 5 6 5 3 6 5 5 5 5 ...
 $ type_accident         : num [1:20065] 2 2 7 2 2 7 3 3 3 2 ...
 $ accident_car          : num [1:20065] 1 1 0 1 1 0 1 1 1 1 ...
 $ accident_passenger    : num [1:20065] 0 0 1 0 0 1 0 0 0 0 ...
 $ accident_motor        : num [1:20065] 0 0 0 0 0 0 0 0 0 0 ...
 $ accident_goodsroad    : num [1:20065] 0 0 0 0 0 0 0 0 0 0 ...
 $ accident_other        : num [1:20065] 0 0 0 0 0 0 0 0 0 0 ...
 $ LINREFX               : num [1:20065] 582353 565114 567468 562842 555294 ...
 $ LINREFY               : num [1:20065] 5926138 5939258 5937258 5941421 5936479 ...
 $ road_surface_condition: num [1:20065] NA NA NA NA NA NA NA NA NA NA ...
 $ day                   : num [1:20065] 25 10 18 5 15 17 18 8 26 5 ...
 $ X_fid                 : chr [1:20065] "DE.HH.UP_RADWEGENETZ_GESAMT_111497" "DE.HH.UP_RADWEGENETZ_GESAMT_12757" "DE.HH.UP_RADWEGENETZ_GESAMT_42600" "DE.HH.UP_RADWEGENETZ_GESAMT_86209" ...
 $ street_name           : chr [1:20065] "Rothenhauschaussee" "Tarpenbekstraße" "Herderstraße" "Garstedter Weg" ...
 $ type_bikelane         : chr [1:20065] "Calle de tráfico mixto > 50 km/h" "Calle de tráfico mixto ≤ 30 km/h" "Otros" "Otros" ...
 $ direction             : chr [1:20065] "in Geometrie-Richtung" "in beide Richtungen" "in Geometrie-Richtung" "in beide Richtungen" ...
 $ surface               : chr [1:20065] "Bituminöse Decke" "Bituminöse Decke" "Kunststein-Pflaster" "Bituminöse Decke" ...
 $ width                 : num [1:20065] 3 3.5 1.5 2.3 NA 1.7 2.25 6 1.5 2 ...
 $ pressure              : num [1:20065] 1003 990 1019 996 1006 ...
 $ rainfall              : num [1:20065] 0 0 0 0.1 0 0 0 0 0 0 ...
 $ sunshine_minutes      : num [1:20065] 0 0 0 0 0 60 0 1 60 0 ...
 $ temperature           : num [1:20065] 3.7 3.9 5.8 4.1 1.2 -0.6 -1.1 5.7 14.1 13.2 ...
 $ humidity              : num [1:20065] 91 90 84 94 94 92 95 77 55 88 ...
 $ visibility            : num [1:20065] 18000 14000 30000 4800 7000 12000 2400 20000 40000 20000 ...
 $ weather_description   : chr [1:20065] "nieve" "lluvia" "desconocido" "lluvia" ...
 $ wind_direction        : chr [1:20065] "SW" "SW" "W" "NE" ...
 $ wind_category         : chr [1:20065] "Fresco" "Fresco" "Fresco" "Muy débil" ...
 $ traffic               : num [1:20065] 138 193 258 69 491 407 245 163 549 361 ...

Además realizaremos un análisis estadísticos de los datos mediante summary(df)

summary(df)
      ...1            FID         admin_region   admin_district  
 Min.   :    1   Min.   : 6188   Min.   :1.000   Min.   : 1.000  
 1st Qu.: 5017   1st Qu.:95211   1st Qu.:2.000   1st Qu.: 3.000  
 Median :10033   Median :96654   Median :3.000   Median : 6.000  
 Mean   :10033   Mean   :87810   Mean   :3.328   Mean   : 6.693  
 3rd Qu.:15049   3rd Qu.:98070   3rd Qu.:5.000   3rd Qu.: 9.000  
 Max.   :20065   Max.   :99950   Max.   :7.000   Max.   :20.000  
                 NA's   :17859                                   
  municipality        year          month             hour      
 Min.   :101.0   Min.   :2016   Min.   : 1.000   Min.   : 0.00  
 1st Qu.:212.0   1st Qu.:2018   1st Qu.: 5.000   1st Qu.: 9.00  
 Median :319.0   Median :2020   Median : 7.000   Median :14.00  
 Mean   :348.2   Mean   :2020   Mean   : 6.898   Mean   :13.33  
 3rd Qu.:507.0   3rd Qu.:2022   3rd Qu.: 9.000   3rd Qu.:17.00  
 Max.   :718.0   Max.   :2023   Max.   :12.000   Max.   :23.00  
                                                                
    day_week        category     kind_accident   type_accident  
 Min.   :1.000   Min.   :1.000   Min.   :0.000   Min.   :1.000  
 1st Qu.:3.000   1st Qu.:3.000   1st Qu.:1.000   1st Qu.:2.000  
 Median :4.000   Median :3.000   Median :4.000   Median :3.000  
 Mean   :4.072   Mean   :2.906   Mean   :3.262   Mean   :4.024  
 3rd Qu.:6.000   3rd Qu.:3.000   3rd Qu.:5.000   3rd Qu.:6.000  
 Max.   :7.000   Max.   :3.000   Max.   :9.000   Max.   :7.000  
                                                                
  accident_car    accident_passenger accident_motor     accident_goodsroad
 Min.   :0.0000   Min.   :0.00000    Min.   :0.000000   Min.   :0.0000    
 1st Qu.:0.0000   1st Qu.:0.00000    1st Qu.:0.000000   1st Qu.:0.0000    
 Median :1.0000   Median :0.00000    Median :0.000000   Median :0.0000    
 Mean   :0.5861   Mean   :0.07506    Mean   :0.009569   Mean   :0.0145    
 3rd Qu.:1.0000   3rd Qu.:0.00000    3rd Qu.:0.000000   3rd Qu.:0.0000    
 Max.   :1.0000   Max.   :1.00000    Max.   :1.000000   Max.   :1.0000    
                                                        NA's   :2180      
 accident_other     LINREFX          LINREFY        road_surface_condition
 Min.   :0.000   Min.   :548947   Min.   :5917209   Min.   :0.000         
 1st Qu.:0.000   1st Qu.:563047   1st Qu.:5934205   1st Qu.:0.000         
 Median :0.000   Median :565977   Median :5936362   Median :0.000         
 Mean   :0.088   Mean   :566542   Mean   :5936301   Mean   :0.269         
 3rd Qu.:0.000   3rd Qu.:569844   3rd Qu.:5939078   3rd Qu.:0.000         
 Max.   :1.000   Max.   :586908   Max.   :5952865   Max.   :2.000         
 NA's   :13356                                      NA's   :11526         
      day           X_fid           street_name        type_bikelane     
 Min.   : 1.00   Length:20065       Length:20065       Length:20065      
 1st Qu.: 8.00   Class :character   Class :character   Class :character  
 Median :16.00   Mode  :character   Mode  :character   Mode  :character  
 Mean   :15.72                                                           
 3rd Qu.:23.00                                                           
 Max.   :31.00                                                           
                                                                         
  direction           surface              width           pressure     
 Length:20065       Length:20065       Min.   : 0.750   Min.   : 977.6  
 Class :character   Class :character   1st Qu.: 1.850   1st Qu.:1010.1  
 Mode  :character   Mode  :character   Median : 3.000   Median :1016.2  
                                       Mean   : 3.352   Mean   :1015.5  
                                       3rd Qu.: 5.000   3rd Qu.:1021.3  
                                       Max.   :16.000   Max.   :1047.9  
                                       NA's   :141      NA's   :16      
    rainfall       sunshine_minutes  temperature        humidity    
 Min.   : 0.0000   Min.   : 0.00    Min.   :-12.70   Min.   : 20.0  
 1st Qu.: 0.0000   1st Qu.: 0.00    1st Qu.:  8.30   1st Qu.: 57.0  
 Median : 0.0000   Median :11.00    Median : 14.60   Median : 72.0  
 Mean   : 0.0931   Mean   :23.06    Mean   : 13.97   Mean   : 70.8  
 3rd Qu.: 0.0000   3rd Qu.:52.00    3rd Qu.: 19.30   3rd Qu.: 86.0  
 Max.   :13.5000   Max.   :60.00    Max.   : 37.70   Max.   :100.0  
 NA's   :71        NA's   :1279     NA's   :14       NA's   :14     
   visibility    weather_description wind_direction     wind_category     
 Min.   :   80   Length:20065        Length:20065       Length:20065      
 1st Qu.:20000   Class :character    Class :character   Class :character  
 Median :39640   Mode  :character    Mode  :character   Mode  :character  
 Mean   :34117                                                            
 3rd Qu.:45000                                                            
 Max.   :80000                                                            
 NA's   :6                                                                
    traffic      
 Min.   :   0.0  
 1st Qu.: 254.0  
 Median : 418.0  
 Mean   : 474.7  
 3rd Qu.: 628.0  
 Max.   :2256.0  
 NA's   :24      

Para realizar una matriz de covarianzas y el análisis de PCA, las variables deben de ser numéricas. Tal y como hemos visto anteriormente, existen algunas variables de tipo carácter, que a continuación transformaremos a tipo factor y finalmente a tipo numérico para poder implementar tanto la matriz de covarianzas como el análisis de PCA.

# Convertimos las columnas en numericas
df$category <- as.factor(df$category)
df$category <- as.numeric(df$category)

df$type_bikelane <- factor(df$type_bikelane)
df$type_bikelane<- as.numeric(df$type_bikelane)

df$direction <- factor(df$direction)
df$direction <- as.numeric(df$direction)

df$surface <- factor(df$surface)
df$surface <- as.numeric(df$surface)

df$weather_description <- factor(df$weather_description)
df$weather_description <- as.numeric(df$weather_description)

df$wind_direction <- factor(df$wind_direction)
df$wind_direction <- as.numeric(df$wind_direction)

df$wind_category <- factor(df$wind_category)
df$wind_category <- as.numeric(df$wind_category)

Descartamos además las siguientes 3 variables que no son de interés para la creación de los histogramas a partir de los cuales se puede ver la distribución de las diferentes variables.

df <- df %>% select(-...1, -FID, -X_fid, -street_name)

2. Análisis de la distribución de las diferentes variables

Para realizar el análisis estadístico de las diferentes variables que describen a los accidentes ciclistas en Hamburgo, en primer lugar crearemos un histograma para cada una de las variables:

for (col in names(df)) {
  
  if (col == "category") next
  
  variable <- df[[col]]
  
  if (is.numeric(variable)) {
    
    df_filtered <- df 
    
    # Agrupamos en bins
    bin_data <- df_filtered %>%
      mutate(bin = cut(!!sym(col), breaks = 31, include.lowest = TRUE)) %>%
      group_by(bin, category) %>%
      summarise(n = n(), .groups = "drop")
    
    # Ordenamos los niveles de los bins
    bin_data$bin <- factor(bin_data$bin, levels = unique(bin_data$bin))
    
    # Creamos el histograma
    p <- ggplot(bin_data, aes(x = bin, y = n, fill = as.factor(category))) +
      geom_bar(stat = "identity") +
      labs(title = paste("Cantidad de accidentes por", col, "y category"),
           x = col, y = "Número de accidentes", fill = "category") +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
    
    print(p)
  }
}

En los histogramas representados se pueden observar las siguientes tendencias dentro de los accidentes ciclistas:

  • Existen ciertas zonas en la ciudad de Hamburgo donde se concentran los accidentes ciclistas.
  • La cantidad de accidentes, especialmente los accidentes leves, han aumentado en los últimos años.
  • La cantidad de accidentes aumenta en los meses cálidos, probablemente debido al mayor tráfico.
  • Durante el día, los accidentes se concentran en la mañana (7-8 horas) cuando la gente va al trabajo y a la hora de salida del trabajo (16-17 horas), probablemente debido al mayor tráfico en estas horas.
  • La cantidad de accidentes disminuye el fin de semana.
  • En la mayoria de accidentes ciclistas también se ve involucrado un coche pero no peatones.
  • A medida que aumenta la humedad, aumenta la cantidad de accidentes.
  • En relación al tráfico, los accidentes aumentan a medida que el tráfico aumenta, pero una vez alcanzado un máximo, la cantidad de accidentes disminuye.

Además realizamos un análisis estadístico con las variables normalizadas en función del tráfico:

for (col in names(df)) {
  
  if (col == "category" || col == "traffic") next
  
  variable <- df[[col]]
  
  if (is.numeric(variable)) {
    
    df_filtered <- df 
    
    # Agrupamos en bins y sumamos tráfico
    bin_data <- df_filtered %>%
      mutate(bin = cut(!!sym(col), breaks = 31, include.lowest = TRUE)) %>%
      group_by(bin, category) %>%
      summarise(
        n = n(),
        total_traffic = sum(traffic, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      mutate(rate = n / total_traffic)  # Normalizamos
    
    # Ordenamos los niveles de los bins
    bin_data$bin <- factor(bin_data$bin, levels = unique(bin_data$bin))
    
    # Creamos el histograma con valores normalizados
    p <- ggplot(bin_data, aes(x = bin, y = rate, fill = as.factor(category))) +
      geom_bar(stat = "identity") +
      labs(title = paste("Tasa de accidentes por", col, "y category (normalizado por tráfico)"),
           x = col, y = "Tasa de accidentes (por unidad de tráfico)", fill = "category") +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
    
    print(p)
  }
}

Los histogramas normalizados indican lo siguiente:

  • La cantidad de accidentes en función del tráfico se mantiene relativamente constante.
  • Los accidentes en función del tráfico aumentan en los meses más cálidos.
  • Existen tipos de carriles donde los accidentes mortales no ocurren o ocurren mucho menos que en otro tipo de carriles.

3. Matriz de correlaciones

En esta sección calcularemos y visualizaremos la matriz de correlaciones entre las variables numéricas del dataset. Esto nos permitirá identificar relaciones lineales entre variables, lo cual nos podrá ser útil para la selección posterior de variables en los modelos predictivos.

Primero seleccionaremos solo las columnas numéricas del dataframe.

numeric_df <- df[sapply(df, is.numeric)]

Para crar la matriz de correlaciones correctamente, deberemos de tener un valor razonable de NAs. Por ello, vamos a controlar primero cuantos elementos NA hay para cada columna.

colSums(is.na(numeric_df))
          admin_region         admin_district           municipality 
                     0                      0                      0 
                  year                  month                   hour 
                     0                      0                      0 
              day_week               category          kind_accident 
                     0                      0                      0 
         type_accident           accident_car     accident_passenger 
                     0                      0                      0 
        accident_motor     accident_goodsroad         accident_other 
                     0                   2180                  13356 
               LINREFX                LINREFY road_surface_condition 
                     0                      0                  11526 
                   day          type_bikelane              direction 
                     0                    141                    141 
               surface                  width               pressure 
                   141                    141                     16 
              rainfall       sunshine_minutes            temperature 
                    71                   1279                     14 
              humidity             visibility    weather_description 
                    14                      6                     77 
        wind_direction          wind_category                traffic 
                    88                    116                     24 

A continuación, definiremos un umbral de proproción mínima de datos correctos y filtraremos las columnas numéricas según este umbral, con tal de descartar las columnas con muchos valores NAs.

# Definimos un umbral para la proporción mínima de datos no faltantes (80%)
threshold <- 0.8

# Filtramos las columnas numéricas que tienen al menos el 80% de valores no nulos
valid_cols <- sapply(numeric_df, function(x) mean(!is.na(x)) > threshold)
numeric_df_filtered <- numeric_df[, valid_cols]

Una vez realizado esto, podemos calcular la matriz de correlación y representarla graficamente:

# Calculamos la matriz de correlación
cor_matrix <- cor(numeric_df_filtered, use = "complete.obs")

# Representamos la matriz
corrplot(cor_matrix, method = "color", type = "upper", 
         tl.col = "black", tl.cex = 0.8, number.cex = 0.7,
         col = colorRampPalette(c("blue", "white", "red"))(200))

En la matriz de correlaciones, vemos que existe una correlación fuerte entre las diferentes variables climatológicas consideradas, al igual que existe una correlación entre las variables espaciales y los limites administrativos.

En cuanto al tipo y la categoria de accidente, se puede apreciar cierta correlación con alguna variable climatológica (humedad, temperatura) al igual que con el tipo de carril ciclista.

Para posteriormente utilizar el dataframe numeric_df en los modelos de aprendizaje automático, a continuación guardamos los datos en un csv.

### Guardamos el numeric_df para posteriormente utilizarlo en los modelos
write.csv(numeric_df, "data_complete_numeric.csv", row.names = FALSE)

4. Análisis de PCA

En esta sección se realiza un Análisis de Componentes Principales (PCA) para reducir la dimensionalidad del conjunto de datos numéricos y facilitar la interpretación visual de las relaciones entre variables. El PCA permite identificar las combinaciones lineales de variables que explican mayor varianza en los datos, lo que puede ser útil tanto para exploración como para selección de variables relevantes.

Antes de aplicar el PCA, se realiza una limpieza del conjunto de datos. Para ello, eliminamos las columnas que contienen valores no finitos (Inf) o valores faltantes (NA), ya que estos afectarían negativamente el análisis.

# Eliminamos las columnas que contienen NA o Inf
numeric_df_clean <- numeric_df_filtered %>%
  select(where(~ all(is.finite(.)) & !any(is.na(.))))

Por otra parte escalamos las variables numéricas para que todas tengan media cero y desviación estándar uno. Esto es fundamental en PCA, ya que de lo contrario, las variables con mayores magnitudes dominarían el análisis, sin necesariamente ser más importantes.

# Escalamos las variables
numeric_df_scaled <- scale(numeric_df_clean)

Una vez realizado estos pasos, aplicamos la función prcomp para calcular el PCA. Además extraemos las cargas (loadings), que indican el peso de cada variable en cada componente principal.

# Realizamos el PCA
pca_result <- prcomp(numeric_df_scaled, center = TRUE, scale. = TRUE)

# Cargas absolutas
loadings <- abs(pca_result$rotation[, 1:3])

A continuación extraemos las cargas absolutas de las tres primeras componentes principales y identificamos las cinco variables más influyentes en cada una de ellas. Estas son las variables que más peso tienen en la definición de cada componente, y por tanto son claves para entender qué representa cada uno.

El resultado top_vars muestra estas cargas, ordenadas por importancia:

# Seleccionamos las 5 variables más influyentes por componente
top_vars <- apply(loadings, 2, function(x) sort(x, decreasing = TRUE)[1:5])

# Mostramos resultados
top_vars
           PC1       PC2       PC3
[1,] 0.6210605 0.5671171 0.7384767
[2,] 0.6170232 0.5350412 0.4721833
[3,] 0.4165970 0.5058274 0.3187590
[4,] 0.1578569 0.3019535 0.1941320
[5,] 0.1083283 0.1053220 0.1775266
# Seleccionamos los nombres de las 5 variables más influyentes por componente
top_var_names <- apply(loadings, 2, function(x) names(sort(x, decreasing = TRUE)[1:5]))

# Mostramos resultados
top_var_names
     PC1              PC2             PC3                 
[1,] "municipality"   "type_accident" "accident_passenger"
[2,] "admin_region"   "accident_car"  "kind_accident"     
[3,] "LINREFX"        "kind_accident" "accident_car"      
[4,] "admin_district" "year"          "LINREFY"           
[5,] "accident_car"   "admin_region"  "admin_district"    
# Seleccionamos los nombres de las 5 variables menos influyentes por componente
bottom_var_names <- apply(loadings, 2, function(x) names(sort(x, decreasing = FALSE)[1:5]))

# Mostramos resultados
bottom_var_names
     PC1              PC2              PC3            
[1,] "day"            "admin_district" "type_accident"
[2,] "month"          "day_week"       "municipality" 
[3,] "day_week"       "day"            "day_week"     
[4,] "accident_motor" "month"          "admin_region" 
[5,] "LINREFY"        "hour"           "month"        
fviz_pca_biplot(pca_result,
                label = "var",       
                habillage = numeric_df_clean$category,
                addEllipses = TRUE,   
                palette = "Dark2",    
                repel = TRUE,         
                col.var = "black",
                title = "Biplot del PCA coloreado por tipo de accidente")

Este gráfico nos muestra tanto:

  • Las observaciones proyectadas sobre los dos primeros componentes principales.
  • Las variables originales representadas como vectores, indicando su dirección y fuerza en el espacio de componentes.

Los puntos han sido coloreados según la variable categórica (category), que en este caso representa el tipo de accidente. Esto permite explorar visualmente que los diferentes tipos de accidente se agrupan de forma natural en el espacio reducido de los primeros componentes de forma similar, puesto que las elipses de cada grupo tienen formas relativamente parecidas y los grupos están poco separados entre sí a excepción quizá del grupo 1.

4. Análisis geoespacial

Además, para finalizar el análisis de los datos incluiremos un análisis geoespacial de los datos para visualizar su distribución en el espacio.

df <- read_csv("data_complete.csv", locale = locale(encoding = "UTF-8"))
df_sf <- st_as_sf(df, coords = c("LINREFX", "LINREFY"), crs = 25832) %>% 
  st_transform(4326)

La primera variable que hemos considerado para ello ha sido la categoria del accidente:

category_colors <- colorFactor(c("blue", "red", "green"), domain = levels(factor(df_sf$category)))

leaflet(df_sf) %>%
  addTiles() %>%
  addCircleMarkers(
    ~st_coordinates(df_sf)[, 1], ~st_coordinates(df_sf)[, 2], 
    color = ~category_colors(category),  
    popup = ~category,
    radius = 3, 
    fillOpacity = 0.7
  ) %>%
  addLegend(
    "bottomright", 
    pal = category_colors,  
    values = ~category,     
    title = "Tipo de accidente"
  )

A continuación también representamos la distribución de los accidentes según el tipo de accidente.

type_accident_colors <- colorFactor(
  palette = c("blue", "red", "green", "yellow", "purple", "orange", "pink"),
  domain = levels(factor(df_sf$type_accident))
)

leaflet(df_sf) %>%
  addTiles() %>%
  addCircleMarkers(
    ~st_coordinates(df_sf)[, 1], ~st_coordinates(df_sf)[, 2], 
    color = ~type_accident_colors(type_accident),  
    popup = ~type_accident,
    radius = 3, 
    fillOpacity = 0.7
  ) %>%
  addLegend(
    "bottomright", 
    pal = type_accident_colors,  
    values = ~type_accident,    
    title = "Tipo de accidente"
  )

Finalmente según el tipo de carril bici donde ha ocurrido el accidente.

type_bikelane_colors <- colorFactor(
  palette = c("blue", "red", "green", "yellow", "purple", "orange", "pink", "brown"),
  domain = levels(factor(df_sf$type_bikelane))
)

leaflet(df_sf) %>%
  addTiles() %>%
  addCircleMarkers(
    ~st_coordinates(df_sf)[, 1], ~st_coordinates(df_sf)[, 2], 
    color = ~type_bikelane_colors(type_bikelane),  
    popup = ~type_bikelane,
    radius = 3, 
    fillOpacity = 0.7
  ) %>%
  addLegend(
    "bottomright", 
    pal = type_bikelane_colors,  
    values = ~type_bikelane,     
    title = "Tipo de carril bici"
  )

Como se puede apreciar de las 3 gráficas, no se puede identificar ninguna concentración espacial de los accidentes en el espacio según las 3 variables seleccionadas.

A continuación, y para finalizar el análisis geospacial de los accidentes ciclistas se representa un heatmap que permite identificar visualmente las zonas con mayor concentración de siniestros en la ciudad, facilitando la detección de posibles puntos críticos para la seguridad vial ciclista.

leaflet(df_sf) %>%
  addTiles() %>%
  addHeatmap(
    lng = ~st_coordinates(df_sf)[, 1],
    lat = ~st_coordinates(df_sf)[, 2],
    intensity = 1,       
    blur = 20,
    max = 0.05,
    radius = 15
  )